home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / HyperCard Files / MacMathPascal / card_7528.txt < prev    next >
Encoding:
Text File  |  1988-09-05  |  4.3 KB  |  221 lines
  1. -- card: 7528 from stack: in
  2. -- bmap block id: 0
  3. -- flags: 0000
  4. -- background id: 2607
  5. -- name: LEAST_SQUARE.PAS
  6.  
  7.  
  8. -- part contents for background part 17
  9. ----- text -----
  10. LEAST_SQUARE[procedure]
  11.  
  12.  
  13. PURPOSE 
  14.  
  15. LEAST_SQUARE calculates the least square fit of data to the curve of the form y = Γêæax^n for n = 0 to NMAX. NMAX is less than the number of data  points.
  16.  
  17.  
  18. TYPE REQUIREMENTS
  19. const
  20. MAXPOWER = the maximum power for the polynomial
  21. MAXSIZE = A MAXIMUM SIZE FOR DATA ARRAY MUST BE GREATER THEN MAXPOWER
  22.  
  23. TYPE
  24. FLOAT = REAL or DOUBLE or EXTENDED;
  25. DATA = ARRAY[1..MAX NUMBER OF POINTS] of Float;
  26. RAY = ARRAY [1..NMAX+1] OF FLOAT;
  27. FIXXED = INTEGER of LONGINT;
  28.  
  29.  
  30. CALLING PROCEDURE
  31.  
  32. var
  33.  
  34. X:DATA;{Array containing the independent variable}
  35. Y:DATA;{Array containing the dependent variable}
  36. NDP:FIXXED;{The number of data points}
  37. COEF:RAY;{The output least square parameter of the above equation} 
  38. POWER:FLOAT;{NMAX in the above equation}
  39.  
  40.  
  41. Define X, Y, NDP, POWER.  Then call
  42.  
  43. LEAST_SQUARE(COEF,X,Y,NDP,POWER);.
  44.  
  45. On return COEF contains the coefients of the polynomial least square fit.
  46.  
  47.  
  48.  
  49.  
  50. EXAMPLE
  51.  
  52. The example program reads in NDP, POWER, X, and Y from the keyboard and prints COEF.
  53.  
  54.  
  55. REFERENCES
  56.  
  57. James, M. L., Smith, G. M., Wolford, J. C.:  Applied Numerical Methods for Digital Computation With Fortran and CSMP, Harper and Row, New York, 1977.
  58.  
  59.  
  60. -- part contents for background part 22
  61. ----- text -----
  62. LEAST_SQUARE.PAS
  63.  
  64. -- part contents for background part 23
  65. ----- text -----
  66. LEAST_SQUARE.PAS
  67.  
  68. -- part contents for background part 26
  69. ----- text -----
  70. 18
  71.  
  72. -- part contents for background part 6
  73. ----- text -----
  74. PROGRAM TEST_LEAST_SQUARE (INPUT, OUTPUT);
  75. CONST
  76. MAXSIZE = 10;
  77. maxpower = 9;
  78. TYPE
  79. FLOAT = extended;
  80. FIXXED = INTEGER;
  81. SQUARE = ARRAY[1..MAXSIZE, 1..MAXSIZE] OF FLOAT;
  82. RAY = ARRAY[1..20] OF FLOAT;
  83. DATA = ARRAY[1..100] OF FLOAT;
  84. VAR
  85. X, Y : DATA;
  86. COEF : RAY;
  87. INPUTS, POWER, NDP : FIXXED;
  88. PROCEDURE see;
  89. VAR
  90. R : Rect;
  91. BEGIN
  92. HideAll;
  93. SetRect(R, 0, 35, 550, 330);
  94. SettextRect(R);
  95. Showtext;
  96. END;
  97. PROCEDURE LEAST_SQUARE (VAR COEF : RAY;
  98. X, Y : DATA;
  99. NDP, POWER : FIXXED);
  100. VAR
  101. I, J, K, NCOEF : FIXXED;
  102. A : SQUARE;
  103. LEFT : RAY;
  104. f : ARRAY[1..maxsize, 0..maxpower] OF float;
  105. PROCEDURE INVERT (VAR MATRIX : SQUARE;
  106. SIZE : FIXXED);
  107. VAR
  108. SWITCH : ARRAY[1..maxsize, 1..2] OF FIXXED;
  109. K, JJ, KP1, I, J, L, KROW, IROW : FIXXED;
  110. PIVOT, TEMP : FLOAT;
  111. BEGIN
  112. FOR K := 1 TO SIZE DO
  113. BEGIN
  114. JJ := K;
  115. IF (K <> SIZE) THEN
  116. BEGIN
  117. KP1 := K + 1;
  118. PIVOT := ABS(MATRIX[K, K]);
  119. FOR I := KP1 TO SIZE DO
  120. BEGIN
  121. TEMP := ABS(MATRIX[I, K]);
  122. IF (PIVOT < TEMP) THEN
  123. BEGIN
  124. PIVOT := TEMP;
  125. JJ := I;
  126. END;
  127. END;
  128. END;
  129. SWITCH[K, 1] := K;
  130. SWITCH[K, 2] := JJ;
  131. IF (JJ <> K) THEN
  132. FOR J := 1 TO SIZE DO
  133. BEGIN
  134. TEMP := MATRIX[JJ, J];
  135. MATRIX[JJ, J] := MATRIX[K, J];
  136. MATRIX[K, J] := TEMP;
  137. END;
  138. FOR J := 1 TO SIZE DO
  139. IF (J <> K) THEN
  140. MATRIX[K, J] := MATRIX[K, J] / MATRIX[K, K];
  141. MATRIX[K, K] := 1.0 / MATRIX[K, K];
  142. FOR I := 1 TO SIZE DO
  143. IF (I <> K) THEN
  144. FOR J := 1 TO SIZE DO
  145. IF (J <> K) THEN
  146. MATRIX[I, J] := MATRIX[I, J] - MATRIX[K, J] * MATRIX[I, K];
  147. FOR I := 1 TO SIZE DO
  148. IF (I <> K) THEN
  149. MATRIX[I, K] := -MATRIX[I, K] * MATRIX[K, K];
  150. END;
  151. FOR L := 1 TO SIZE DO
  152. BEGIN
  153. K := SIZE - L + 1;
  154. KROW := SWITCH[K, 1];
  155. IROW := SWITCH[K, 2];
  156. IF (KROW <> IROW) THEN
  157. FOR I := 1 TO SIZE DO
  158. BEGIN
  159. TEMP := MATRIX[I, KROW];
  160. MATRIX[I, KROW] := MATRIX[I, IROW];
  161. MATRIX[I, IROW] := TEMP;
  162. END;
  163. END;
  164. END;
  165. BEGIN
  166. {define f}
  167. FOR k := 1 TO ndp DO
  168. BEGIN
  169. f[k, 0] := 1.0;
  170. IF (x[k] <> 0.0) THEN
  171. BEGIN
  172. FOR i := 1 TO power DO
  173. f[k, i] := f[k, i - 1] * x[k];
  174. END
  175. ELSE
  176. BEGIN
  177. FOR i := 1 TO power DO
  178. f[k, i] := 0.0;
  179. END;
  180. END;
  181. NCOEF := POWER + 1;
  182. FOR I := 1 TO NCOEF DO
  183. FOR J := 1 TO NCOEF DO
  184. BEGIN
  185. A[I, J] := 0.0;
  186. FOR K := 1 TO NDP DO
  187. A[I, J] := A[I, J] + F[K, I - 1] * F[K, J - 1];
  188. END;
  189. INVERT(A, NCOEF);
  190. FOR I := 1 TO NCOEF DO
  191. BEGIN
  192. LEFT[I] := 0.0;
  193. FOR K := 1 TO NDP DO
  194. LEFT[I] := LEFT[I] + F[K, I - 1] * Y[K];
  195. END;
  196. FOR I := 1 TO NCOEF DO
  197. BEGIN
  198. COEF[I] := 0.0;
  199. FOR K := 1 TO NCOEF DO
  200. COEF[I] := COEF[I] + A[I, K] * LEFT[K];
  201. END;
  202. END;
  203. BEGIN
  204. see;
  205. WRITELN(' HOW MANY DATA POINTS ?');
  206. READLN(NDP);
  207. WRITELN(' WHAT ORDER OF FIT ? ( MUST BE LESS THEN ', NDP, ')');
  208. READLN(POWER);
  209. FOR INPUTS := 1 TO NDP DO
  210. BEGIN
  211. WRITELN(' ENTER X AND Y');
  212. READLN(X[INPUTS], Y[INPUTS]);
  213. END;
  214. LEAST_SQUARE(COEF, X, Y, NDP, POWER);
  215. WRITELN('THE LEAST SQUARE FIT IS - Y =');
  216. FOR INPUTS := 1 TO (POWER + 1) DO
  217. IF (INPUTS = 1) THEN
  218. WRITELN(' ', COEF[INPUTS] : 20 : 10)
  219. ELSE
  220. WRITELN('+', COEF[INPUTS] : 20 : 10, ' * X  TO THE ', (INPUTS - 1) : 3);
  221. END.